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ABSTRACT 

With the advent of the high speed digital computer, many 
problems heretofore considered unsolvable for all practical 
purposes are now well within the reach of the anplied math- 
ematician. One such problem is the routing of a shin through 
a time dependent ocean wave field, from one point on the 
earth's surface to another, so as to minimize a cost func- 
tion of the storm g(x,y, t, a). 

This paper considers a numerical solution to the above 
problem. The technique to be employed is known as the method 
of steepest ascent and is attributed to Arthur E. Bryson and 
Walter F. Denham [1]. Although the computer program as given 
in the Appendix is written specifically for a VC2AP3 class 
vessel operating in a described area of the North Pacific 
Ocean, it can be readily modified to accommodate any type 


vessel operating in the Northern Hemisphere. 
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a, Statement of the Problem 

The basic problem to be developed and solved in this 
paper is the following: A VC2AP3 vessel, located at latitude 
40°N longitude 154°E, is to transit a time dependent ocean 
wave height and wave direction field to latitude 38°N longi- 
tude 123°W so as to minimize a cost function of the form 
e6e,y,u,t) . 

A south-polar stereographic projection of the Northern 
Hemisphere upon a plane passing through the circle of 60°N 
latitude is used to establish a rectangular coordinate sys- 
tem, with the OX and OY axes parallel to the projections of 
the meridians of 10°E and 100°E longitude, respectively. 

Then a 63 x 63 grid is constructed along the OX and OY axes, 
with x = y = 31 defining the projection of the North Pole. 

At each point of this grid, actual wave heights and wave 
directions taken from the files of the U.S. Navy Fleet Numer- 
ical Weather Central, Monterey, California, are recorded. 

The map scale factor MSF, defined as the ratio of a differen- 
tial distance in the OXY plane to the corresponding differen- 


tial distance on the earth's surface, is 


MSR =" [973.75+(x-31) 7+ ¢y-31)°1/1043.6. 

Germane to the solution of the above problem is the 
elliptical polar velocity diagram of Figure l, giving the 
ship's speed v as a function of the angle 0 between the 
ship's heading and the wave direction. R. W. James [6] gives 
empirical curves for the speed of a VC2AP3 class vessel in 


head waves Vn? beam waves Vp! and following waves Ve as a 


function of wave height. Bleick [2] showed through a least 
Squares analysis that these three speeds can be represented 
closely by the 4-parameter formula 


C 


— 2 — 
ve C, (HtC, ) exn ( CH) . 


where H is the wave height. 

Thus, the projected speed of the ship at any point on 
the grid is 

Vivo) t) = v(h,O) MSE | 

where v(H,0) is the ship's speed of the polar velocity dia- 
gram of Figure 1, and H(x,y,t) is the intervolated value of 
the Fleet Numerical Weather Central wave height data. If 
we define K(x,y,t) to be the stereographic grid wave direc- 
tion measured counterclockwise from the OX axis, and u(t) 
to be the ship's heading as measured counterclockwise from 
the OX axis, then 


[v sin(u-K) /b]* + {{ct+v cos (u-K \1/a}? or 


ee mee 


b(a® -C cm 0+ (ar =e 2) sin 20] 


v(H,0) = ————. / (Pe cosO + is 2cos* 


where 0 = u-K, a and b are the semi-principal axes of Figure 
1, and c is the distance to the eccentric pole o. We note 
at this point that an elliptical polar velocity diagram, 
such as Figure 1, must be specified for each wave height, 
and that a, b, c are functions of H(x,v,t). Thus, we can 
write the equations of motion of the ship's projection in 
the OXY plane as: 


x(t) WO Sault t) 


II 


y (t) 


SS aaa) 


weibetq ARTOOTeA AeTO, 


°[T eanbty 





eee 


where the (*) over a variable indicates its time derivative, 
and V = a = V(x,y,u,t). Hence, our problem is to choose 
the control angle, u(t), at each point of the transit so as 


to minimize the cost function g. 








II. Abstract of the Theory 

Consider a system of differential equations of the form 
i, a ai ai 
= =t ol»... neewherestte= £u(X,Ue tb) are, functions sof 
class ae X is an n-dimensional vector, U is an m-dimensional 
vector, and t, the independent variable, is a scalar. We 
desire the solution to the above system of equations which 
is optimal in the sense of maximizing (minimizing) the ter- 


minal value of some cost function g(X,t) while simulta- 


cH 
neously satisfying terminal constraints, 


h,(X,t),_7 = 0 (i=l,...N; 0 < N <n). 


The solution of a problem of the above form by the method of 
steepest ascent is concerned with the development of the m 
control variables, uy (t),...u,(t). Essentially the problem 
is solved in two parts. In part I we are interested in 
attaining admissibility, that is satisfying the terminal 
constraints 


h, (X,t) = OG NGO, < N <%n). 


t=T 
In part II we are concerned with the problem of maximizing 
(minimizing) the terminal value of the cost function g. The 
procedure is as follows: We initially guess nominal values 
ef the control, U: in general, the resulting solution will 
neither be admissible nor an extremal. Changes in the con- 
trol are then generated which oe drive us to admissibil- 
ity, that is, U(t) is replaced by U(t) + 6U(t). This pro- 
cedure is continued until we attain admissibility. The ad- 


missible solution, in general, will not be an extremal; so 


we once again generate changes in the control which tend to 

drive us toward an extremal without affecting admissibility. 
The author realizes the above is rather abbreviated; 

however, this discussion is intended to give a brief idea of 


Wheat 1S to LoOlLow. 


U6, 








IIIf. Development of the Theory 
Notation 
Let us start by considering the following system of 


ordinary nonlinear differential equations: 


CS =x = V cos u = ft (x) x7 ,u,t) 
see) x = y = "Vestn u = Be ise Oe 
ve Z = g(x,y,t) = ee ee 


where V = V(x,y,u,t), u = u(t), and g is the cost function 
to be minimized. We can write this as 

xt = etcxl (ty), x7 (t), ult), t), (i = 1,2,3), 0<t <7, 
where t, the independent variable, is a scalar. It will be 
convenient in the development to follow to represent these 


sets of variables by matrices or by vectors, so let us 


define: 
x(t) rere 
X = x? (t) ; F = Coc oc ; = [ue] 
2 (t) 22 (RO yt 
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where U is a single-rowed column matrix since there is but 
one control variable. 

Equations (3.1) then become 

(m2) ie eee 
The vector X is Cablcderac Stare Vecror, Or Stace, amd the 
one dimensional vector U, GhemcomesneL vector aor control. 


We will assume the control to be a continuous function of 


time. 


Lil 


Now the problem as posed in section I can be interpreted 
as follows: We desire to find a curve C, a solution to (3.2) 
for some choice U, which is optimal in the sense of minimizing 
the terminal value of the cost function g while simultaneously 
satisfying the terminal constraints 


Ll 2 


he 
er & Ciayee= 


x (T) =x eo! 
where x and i are the final values of the state correspond- 
ing to the rectangular coordinates of the terminal point in 
the OXY plane. 
The Adjoint System of Equations 

Of primary interest is the relationship between varia- 
tions, Su(t), im the control at any time, t (0 <"t < 1) 
the resulting terminal variations of the state variable, 


SX(T). If we consider (3.2) in component form, the varia- 


tional relationships for the state variables can be expressed 


as: 
ils ae IL iL 
Sx = of bx + ae spas 6 2. Sx> + of. 6u 
aX ox ox au 
2 2 2 2 
ame es ge ee ee ee 
dk 2 3 
aX ox ox du 
3 3 c) 3 
8° = of bx + of 8x? + of &x> + of é6u 
iL 2 3 
ox ox OX du 


Writing this in matrix form we have: 


Bs 
ox = GéX + HOU 
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ron 
ere sft gel se of 
du 





yxt ox? ax | 1 
4 SX 
| 2 2 2 2 ; ; 
gu. } of of sa q - | Sf a ee 
Il 2 
ax aX ax du 
sx 
ae 4 
jG ofS 9£> “hoa 
yx ara Agee du | 
| a L Mt 
> 
and ou = Pour} 


Note in the special case where the control U consists of a 
Single variable, u, $U can be considered a single-rowed 
Soren VeCr@r. 


If we now introduce a new set of variables 


ate 


Lagrange multipliers, which are undefined functions of t, and 
multiply (3.4) through by ae where at is the transvose of 

A, we obtain the scalar relationship 

(3.5) NSS Seeich ee NCIS 


Let us now integrate (3.5) over the interval (0,T) to obtain 


dk 
0 


el 


(3.6) ASX dt = fo AlGéX at + So ATust dt . 


Integrating the integral on the left by parts yields 
T.3)T a pe ee gees eee Tae 
NSC) eee Creche | viGex ty tal, HOU dt 
Rewriting this as 


AY 8X|5 - [G (A7 + A°G) 6k at = fF ATHsU at 


3 





° | 
and choosing A as a solution to he + Ata) 6x = 0 defines the | 
system of equations adjoint to the variational equations (3.4). 


That is, the adjoint system of equations is 


(290) AW = -A'G 
With this Choicesof wh, (3.6) reduces to 

(3.8) AT 8K _ = ATOX|, + [5 ATHSU at 
and, since we are assuming §X(0) = 0, we have a relation be- 
tween the terminal variations of the state variable, éX(T), 
and the variations in the control, SU. Now if A is chosen 
as the specific solution to the adjoint system of equations, 


(3.7)> such that atetime t = T, A(T) = 990,00) ern Geers 


becomes 

6x" (T) = [7 ATHED at 
giving the terminal variation of x" due to variations sU(t). 
Similarly, choosing A, (T) such that the ie? component is l 


and all other components are 0 yields 


(3.9) sxtcpy = fT at 


"HOU at 
2 


We may consider the matrix product \vH as a vector, say A(t), 





when there is more than one control varrable. Then (3.99 Gag 


be written as 


ey ote ede 


and for a maximum change in x (T) we would choose $U (t) par= 
allel to Pe, for 0 <t< T. Courant [3] page’ 223, remem 


= al eS ee 
LOVEme VECEOr A, (t) as the gradient or x, of A. = Vx. 
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Special Variations 
From the previous discussion we have seen for a maximum 
change in xt we would like 6U(t) to be parallel to Fat Or 
Dies & < 1. Seiims motivates the following choice of a swnecial 


Warlattoen Giuche control, 


beSieds0.) SU(t) = e, Vx" (t) fa e Vx" (t) 


It might be worth noting we do not consider x? since at this 
juncture we are only interested in satisfying the terminal 
constraints 


Ll lL 


page oy ae Ker x? (t) as 


emipstituting (3.10) tnto (3.9), and recalling \vH = Vx", we 


obtain 
Ik Ms Tel oe i 1 a 
ope ly e) fi Vx Vx" dt + e, ho Vx XK aaelit 
oo -11) 
2 it Lee AE 2 2 
6x (t) = ey lo Vx ¢ Vx" dt + e IG Vx + Vx* at 


We now define 


so that (3.11) becomes 
; 2 oye 
Sxa(D)e— a) e.7 Git 2) 
For small perturbations, the value of T will be changed 
only a small amount aT, so that 


axt(t) = éxt(r) + xt(r) aT (i=1,2) 


5 


We can express this relationship in matrix form as 


r 7 7 | 
Ama zi xt ey | Ke - xt (pT) 
(iD) | = | 
alamo 2 | » oa Ones 
i ab. | Se £ 


or ZE = AX. Let us now digress for a moment, and recall that 
if we have an overdetermined system of equations, say AX = B, 
where A, the coefficient matrix, is 3 x 2, X, the unknown 
matrix, 1S 2 x 1, and B is a 3 x l matrix, then a least 
squares solution for X is obtained by multiplying both sides 
of the above equation on the left by Nes provided AA’ 1s net 
Singular. Since we have a system of two equations in three 


unknowns, let us apply the above technique in reverse. 


C 
il 
Let E = zéc where C = so that (3.12) 
C 
2 
becomes zz'c = AX. Now Zz) is a 2 X 2 Gram matrix, so let 
oo ae 
Lia = a= : 
oe 7 


Thus (3.12) is representable as AC = AX, and we are now able 


to solve for C. Use of (3.10) and E = atc gives us our var- 





Lac lLons  1mMeconeRo | sU(t), Reis 10 mc La 

At this point in the development we should be able to 
obtain an admissible curve. We will now consider the pvrob- 
lem of decreasing the cost function without affecting x (D) 
and x? (T). To accomplish this end we consider a special 


variation of the form 


UG 


(3.13) 6U(t), = ay Fx* (t) a ri youn Ge 


2 5 


Sulestitutings (3.13) into (3.9) yields 


3 re 
) 044 j and 
j=l rey.) 


sxt(t) + xt(t) at 


$x2 (TP) 


dxt (T) 


In matrix notation this becomes 


a a) Ee ae + 


f “a 
(3.14) | ; 
2 2H) gel 2m ap? eg ke eh, ; 
f | 2 
4 Sak 5 all MRCS a 
| 
aT 
™ 


Since we have an admissible path, 
ce I xe eRe ep ~ ON 


3 3 e e e 
however, Xe 7 x (T) is unknown, and, in general, is not zero. 


Once again we are faced with the problem of solving an under- 
determined system of equations, (3.14). As before we will 


let a = zip where D = | so that (3.14 becomes AX = Ze. 





Now ZZ) is a 3 x 3 symmetric Gram matrix so let 











Thus (3.14) is representable as 


(3.15) ae =2eD 


Define d. = EB. 3 (1=1,2,3), where € is to be determined, 


and Bs 3 as the cofactor of b. 3 in the determinant of B, that 


is 
d, = elby, bo3 - Boy 443] 
do = €[by, 43 - by, %93! 
d, = €[by, Pop - Pyo Py2] 


Then, upon substieutioen into... Loy, we -omcalit 


xe ~ xt (Tp) 0 
(3.06) oe x4 - x? (T) = "8 0 ‘ 
x = x3 (7) det B 


Thus, by choosing € negative, we can decrease x? and keep xt 


and x? fixed so long as B is nonsingular. That is 


a. “Sees CORand {-— =sdeeeb- 


Now, having D and a = z'p, (3.13) may be solved for SU as a 


multiple of ¢€. In order to limit the maximum change in the 


control, we define the following norm: 


N = {max |u(t)]| + 1 ae! : 
t 
We now define ¢€ = a where DELUMX is the maximum 


allowable change in the control. Thus s(t), the solution 


to (3.13) is replaced by 


€ $U(t) - 
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IV. The Numerical Solution 

The numerical solution of the problem as posed in Sec- 
tion I is based upon the solution of the original system of 
equations, the adjoint equations, and the variational equa- 
tions. For the sake of simplicity, a finite difference 
approximation to the derivative with a variable time step of 
integration is used to integrate the above equations. The 
variable time step of integration is obtained by dividing 
the total time as determined after each path into two hun- 
dred equal increments. Initially, the value for T, the 
terminal time, is determined by the geodesic track, but this 
is not a requirement. Any reasonable guess for T will suf- 
fice. Justification for this method of integration is based 
upon a comparison of the geodesic route as obtained, using 
the above mentioned technique and a fourth order Runge-Kutta 
integration method. 

Due to the core storage limitation imposed by the Fleet 
Numerical Weather Central, 165,000 octal words, a theatre of 
operation is established which is a subfield of the 63 x 63 
grid mentioned in Section I. Having established this sub- 
field the wave field data is then read into core memory. We 
are now in a position to proceed with the numerical solution 
De whe problem. 

Although not a requirement, the geodesic track from the 
point of departure to the point of destination is calculated 


for comparison purposes. As stated in Section II, a nominal 


i? 


control program is initially required. This is taken as the 
angle of inclination of the straight line connecting the 
track end points. 

Recall that in the development of the theory we chose 
as particular solutions to the adjoint equation, Ay: Aor and 


Az, where A, (T) has its Aen component equal to 1 and all other 


components equal to 0. To satisfy this requirement we define 


} 


r 


. 
A, (0) =] 0], A,(0) =; 1], A,(0) =] 0 

0 07 | 1 
and integrate the adjoint system of equations forward to 


obtain A,(T), A(T), A3(T). The resulting matrix 
h, (7) ho (T) A) 
uy (TY yD) 4 (7) 
0,(T) o5(T) o3(T) 


is then inverted to obtain new initial values A, (0), A(0), 


A (0) so that upon integration forward a second time one 


obtains 1 0 | oO 
A, (T) =a Oa \.(T) - ‘al a and A, (T) 7 oe , 
0" co Il! 


Having obtained the initial values of the adjoint vectors 


\ Nos and A. such that at the terminal time T they take on 


dig’ 
prescribed values, we will run a path, check for admissibility, 
and if not admissible generate corrections to the control. 
Simultaneous integration of the original system of equations 


and the adjoint equations gives a path. At each point along 


this path the gradients in the x+ direction (ied, 23) gane 
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calculatedsand stored. for hater neeall.in determining the 
changes in the control. In addition, the elements of the 2Z 
matrix of Section III are calculated. A check is now made 
to see if we have attained admissibility. If not, we pro- 
ceed to calculate the changes in control which will tend to 
satisfy the admissibility constraint. Assuming we do not 
have an admissible path, the A, C, and E matrices of Section 
III are calculated. Thus we obtain new values for the con- 
trol and the terminal time T, that is, U(t) is replaced by 
U(t) + 6U(t) and T is replaced by T + dT. The above proce- 
dure is continued until the admissibility constraint is met. 
Having admissibility we consider the problem of de- 


creasing the cost function 

g = i. e’ Hat ; 
where y is a Scalar, without affecting admissibility. A 
check is made to see if the terminal value of the cost func- 
tion as determined by the last admissible path is less than 
that determined by the previous admissible path. If it is, 
we proceed to calculate the changes in control as defined in 
Section III, so as to decrease the cost function g. If not, 
we decrease the maximum allowable change in the control at 
any time t by one-half. A check is then made to determine 
if the maximum allowable change in control is less than 0.01. 
If it is, a solution has been obtained. If not, proceed to 
calculate the changes in control so as to decrease the cost 
function. Assuming we have no solution we calculate the B, 


D, and a matrices of Section III. This enables one to obtain 


Zag 





new values of the control, that is, U(t) is replaced by 
U(t) + $U(t). A new path is run which, in general, is not 
admissible; so we again drive to admissibility. This proce- 
dure is continued until a solution is obtained. 

It should be noted the values used in the program for 
the maximum allowable change in the control and the stopping 


Criterion are quite arbitrary. 
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Ve Conclusions 

The value of this study lies not in the solution of the 
particular problem considered, but in the fact that given 
any realistic cost function which is a function of the 
stereographic grid coordinates x and y, the control U, and 
time, one can obtain an optimum solution with only slight 
modifications. Such a cost function for a cargo type vessel 
might reflect costs due to storm damage, failure to meet 
scheduled commitment dates, and spoilage of cargo due to ex- 
tended time at sea, to mention but a few. 

On the basis of the results obtained the following 
observations seem justified: 

(a) It was shown in Section III that a neces- 

Sary condition for the solution curve C to be an ex- 

tremal is that the determinant of the B matrix at 

t = T approach zero as track iteration continues. 

In the case where we chose y = 0, that is, minimum 

time, we saw that this, in fact, was the case. In 

the second problem, where we chose y = 0.01, the 

determinant of the B matrix was reduced from a value 

on the order of 5,000 to 0.9. The conclusion to be 

drawn here is that, although the solution curve C 

is not an extremal, it very closely approximates one. 

(b) Defining the initial control as the con- 
SEane LUNnCtT Onsen) = VV (0"°< t <0r), where V re the 
initial angle of departure, leads to convergence to 


within fifteen miles of the point of destination. 
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Thus there seems to be no problem as to how to 
initially choose the control so as to obtain con- 
vergence. 

(c) The simplified integration routine used 
is feasible and yields satisfactory results. 

(d) The requirement of storing an entire con- 
trol program and a complete set of functions 
NH i = 1,2,3 for any one iteration is a marked 
disadvantage. Due to the large storage reaquire- 
ments imposed by the establishment of the wave 
fields a relatively small amount of core memory 
remains. 

(e) An analysis of the results obtained show 
a substantial decrease in the wave heights for the 
minimum time route, as compared to the geodesic 
track, with a decrease in the total transit time 
of approximately 2% hours. The weighted time 
route, on the other hand, shows only slight de- 
creases in the wave heights encountered, as com- 
pared to the minimal time track. This leads one 
to conclude that if the objective is to limit the 
wave heights encountered, a function other than 
exp YH should be used. This is more readily seen 
if we expand at that is, 


em) 


H 
eY = 1+ yH + +5 31 








For small values of y and nominal wave heights, H, 


the dominant factor in minimizing z = Wie eYF at is 


time. 24 
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PROGRAM TEST (INPUT eOUTPIITs TAPES) 

HEMENSTOM TG (12) (200) MEL (200) »GRADX (200) ,GRANY (700), 
1 TES(201) oH4T( 201) 9AZ7(201) 9 AL (201) eWA(201) » WE (39HK9) oWn (63463) 
re XtAM(3) e TLAM(3) 9» XMU(3) 9 XSITGMA(3) eGRADZ (200) ¢ TMII(2) 

COMMON YeVYoebgReCCeHe CK e SK oK x oeKY 9 KT pF IOXeFOYSKOON 

COMMON/L ]ZXAT (R449) « CSK (8448) » SNK (8448) oT 

COMMON/L AJAX HY 

COMMON/SLGESCAPV ¢«CAPVX sCAPVY ee CAP VU sw 

SC OMMON/LAZ XK OXLGe XLT 

PIIITVALFNCE (WE (])ewl(Eel) ) 

wf AL MSF 

ItTEGER NTS 
~ wk witl FIRST ESTABLISH THF wAVF FIELD IN OUR THEATFR OF INTEREST 
TTALXSO F€ NIGX=s0 F IMINXSQN FF JMINKSO F KTX=E0 
' PE i) SXTMENSTONS AN MINT MEIM TNOICES OF ARRAYS 

PEALE GekKX¥eKY eK Te TMIMe IMIMNe CRIT 
Bee ORMAT ( STeeF4,0 ) 
POTNT JekXeKVoKTe IMT Ne IMIN 9 CPIT 
VP PARMAT(JAIHODTMENSTON OF WAVE FTELT ARRAYSS (TL 3elHs T3991 HeT3e) HI) 4X5HT 
Lm INST 34x54 UM INST 3.4 x5HCRI TSF 4,0) 
JOYS IPH Pay"! 
SVS JP em yy 
mee = | ()X 
POy=s JOY 
RCONSK XK Y +k ¥ 
Y¥'FANSFL OATF (KX41) 7/27, 
YME ANSEL OUTE (KY47)77. 
bh YKXEXMP ANN? AN) 
PYysyMF ANP LO] 
ClhF Ay HTS aReay 
rt) ws [TS])ek ft 
Meee 'G({ ys js 
hE ss: OCTAL DATFEFOCTIME GRolpec OF THE TIME SFRIFS MEMRERS 
MER JPLe OTC) eT Tay 
J) fA MAT (COPD eH X ge TA) 
TE CUTGX eEOCITB( YT) cane TTAUX oF OCT TAU ec AND ee TMINKX FO, IMIMe AND, 
PIMINX CEN IMIM CAM KTM OP OKT) 7500675 
GET arkrPaY ATG TO SELFCTEN TIVk SFeETES 
Aes, J! tsk Te j 
me? Lt=leJJ 
Mell TIGR (PTECT) sTTAUeCDTG(T4}) ) 
i Com TTtrue 
Pech SS THE THHIPICKEND aekays 
Pew iN() 4 
Peernn k=l] ,kT 
KKK=( ( Kel )#KYe] ) #KX 
NY 650 KKeJ oe? 
NGsy 
KwokKK=] 
(SLE FIFI ISX (WE eNTG(K) eK Wa NG) 
LF (NGaH1) 90 e990 
44 ChINT G1eMTG(K) 
9S) FORMAT(IGHOFTELN OF DTG=O20e¢) X31THNOT FOUND ON INPUT TAPE4, ARNRT/S) 
GO 10 843 
99 GO TON (97994) KK 
97 [t) & Tal eKx 
LITSKKK+ J 
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jrict=l=hax 
DO &@ JeleKY 
JISJ=KX# Je] I 
DFLYsJUeaJAY 
RKANTSSQORTF ( DELX#NFLX + DELY#NELY ) 
ARG=sWwD (Te TMIN gg J+ JIMIN) /5,72957795) 
CYS=COSF (ARG) 
SINSSTNF (ARG) 
CSK (JJ) = (@=DELX#COS=NEL Y#FSIN) ZROOT 
A SNK (JJ) = (NEL X#STN=DFLY#CNS) JROOT 
GA TH 650 
94 150 J Ts l_kx 
TT= KKK + J] 
iy 19 J=1,KY 
JJ = KXFY + JT 
10 XHT (JJ) Sew CTeIMIN, J@UMIN) 
A501) CONTINUE 
TF CSENSF SwITCRH 1)4109700 
Al POINT J1900¢ KeDTG(K) 
L100 FORMAT (VHONeT2ePH. 9 MPN) 
WO 620 JeleKY 
A620 PPINT 12004 ( wh T*+IMINeJ*UMIN ) og TS19KX ) 
1290 FORMAT( Lxe2hF5S.1/ ) 
TOUS CONT TNE 
PEwIN) 4 
OTGXSNTIG (1) & [LTAUXSTTAU FB TMINXSITMIN &$ JMINXSUMIN $& KTXSKT 
IMO 750 Tal eKkKX 
YO 750 JsleKY 
Sieh kee T 
RISloImMIm % RYJtJeUMIMN F KEREO 
CALL SEALAND (RT,RU LAND, KER) 
JF (RER) 646 4648,646 
A4H PRINT 647 $G$ GN TO 83 
hal FORMAT (26HORT* RI INDEX ERPORS, ARNRT/) 
Aa TF (LAND) 649e750,649 
A49 DON 749 KeleKT 
KKK=((Ke])#KYm])#KX + JU 
749 XHT (KKK) 220.0 
TSG “CONT NUE 
coweF SHOULD NOW HAVE OUR WAVE FIELD DATA STOREN IN MEMMORY 
cr RF AI) HOUR OF NEPARTUPF AND COORDINATES OF TRACK END POINTS 
PFAD 13eHRoeXLG) os XLT1L»9XLGAeXLT2 
13 EEMAT (EF yeily O's ligt De ligt Oe ler oie 1) 
PRINT 14eHPeXLGLoXLT1LoXLG2eXLTe 
14 FORMAT(2]H POUTE OF SHIP BEGINSF4.0,3]H HOURS AFTFR TIME SFRIFS OR 
}IGIN/19H FROM LONGITUDE = F6.191]6H AND LATITIIDE = FAe19/)9H TO 
2 LONGITUDE = F6.et0e16H AND LATITUDF = Fho}) 
roWwF STEEL NOW CONVERT COORDINATES OF TRACK END POINTS Ton GRID VALUES 
beYG=(XIG6)=-10.)/57.42957795) 
COSLG1L=COSF (ARG) 
STNLGLSESINF (ARG) 
LPGS(X| GAa1l,)/57,2957795)] 
COSLG2@=COSF (ARG) 
SINLG2=STNF (ARG) 
ALGEXLT1L/572.2957795) 
COSLTI=ECOSF (ARG) 
SINLTISESINE (ARG) 
APG=XLT?2A/57.27957795) 


28 


mee Ci 


34 


ee 


CASLT2=COSF CARG) 

STNL TASSTNF (ARG) 

FHP S3]P ePNSHCOSLTIZ(1-¢¢SINLT]) 
X}=PRI#COSLG) 

YIL=SPRIFSINLGI 

PRPHZ9)  APNSH#COSL TAS (1. ¢SINLT2) 
XPSPREFCOSLGC 

VP=SPRAHS INI G? 

LSTMAUTayv1l4F ly 

YSTARTHY) oF OY 

ERD SX Pek OX 

YEND2ryrot OY 

UTE GEONESTC rOUTE 

FIESINL TP#COSL TT #STNE GY = COSLTC#HSTINI TI FSTNL GA 
EM= = STmt.Teecose TL¥cOSLG] +* COSLTA*#SINLT1IFCOS] G2 
e‘N= (SINE GOFCOSLG] - COSLGA#SINLG1] ) #COSLT1IHCOSLT2 
PAOTSSGRIF (CELE «¢ FMUEFM + ENZEN ) 
FILSEL/SROOT 

Eee MSE ONT 

ChHeENSROAOT 

MeL XSxPerl 

(FF, Y¥sV¥eey ) 
SUASSARTE (OF LXHDEL X + DEL YOR Y) 
1VC=zS}? 

JF (XLG@ ce Xl) ) 20,27) 92 
AMGSABSF (F N/A? 24) ) 

NNCSASTIUF (ARG#HSI2) /ARG 

r= XS TART 

Y=YSTART 

TSHR/O4, 

»TRES=#0,9 

CALL TERE 

aT ()) sb: 

TFS(CL)SXTFES 

CALL ANGLE 

A7(1LYSEXLG 

pac 1) = Nite T 

oom 1) =X 

2 eae | 

XXSXSTAMT 

YY=SYST4PT 

SS=0.20 

| R= 

STFP=1./Ae 

40 Ke2eP0h 

CIUSP=(FOY"Y) FFN/S3}.2205 * EM 
STNPS(XKeFOX) FENSJ31.2NS |= EL 
CALL TERR 

moet ShTP 

OmblL VFR] v 

X=XX & CAPVECOSPHESTE PB 

YaYY + ChpvVvFtsINpesreEp 

S=SS + CAPVESTEP 

IF (S#AKRC) 359 34434 
KFOTS(AKC.SS)/s(S=@SS) 
TZRATHSTEP + T 

XTEFSSRATHYSTFO + XTES 

Kz (XXX) PPAT «+ XxX 
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)5 


ee 


14 
CWE 


YE(YmYY)#PAT + YY 

Gort One 

=X 

YY =7 

SS2S 

T=T + STEFF 

XTFS=XIFS * STEP 

IF ( (Kel) #2/6 + 1 = NK) 16616615 
VK SNK «@ ] 

C4LL TERP 

MT (NK) oH 

TES (NK) SXTES 

CALL ANGLE 

AZ (NK) Sx¥LG 

AI (NK) =XLT 

wa (NK) =XK 

TF (SwARC) 1494164) 

[F (ARSE (X*XMFAN) PNK) 3795249574 
TF (ARSF (ya yMEAN) eROY) 4095249524 
CTO Ea 

MRINT 7 

FORMAT (////30X14HGENNESTC ROUTE) 
RINT 18 

FORMAT( SX&HDAYSAXSHLONGT4XSHI ATT ASXQGHWAVESXI4HWAVE NIRECTIONs2X9H 
1 OF TRAVELS XSH@ TUDE 4 X4GHTUDES XQAHHE TIGHT 6X) OHF ROM NORTH) 
LQINT Igoe (TES(K) epZ(K) eo AL (K) @HT(K) oWin (K) Kel gNK) 
em AT of moc eP IY loe eR ber Leer is, U ) 
WELL NOW GUFSS AN INITIAL VALUE OF THE CONTROL U(T) AND RUN A PATH, 
DEL TaX=XFNDXSTART 

OF LTAYSYFNDe¥YSTART 

fe CORIETAY «Es. G0.) GOnT ORS 30 
ARGSNEL TAY/SDEL TAX 

ht FASATANF (ARG) 

IF (DELTAX) 53195309932 

TF (DEL TAY) 533,534,534 

VE6,7K319 + ALFA 

GQ 100535 

ve Alor A 

COrr@a 35 

VS3~0e14159 + ALFA 

GO TH $35 

IF (VELTAy)534,537,536 

VE4~e7)124 

GO TO 535 

v2l.570R 

OMe) Sel 

PRINT 539 

FORMAT (Syx4HSTOP) 

G9 TO AR 

CONTINUE 

riIN=Q 

GAMASN,N 

leah? sed. 

WK XSFLOOTF (KX) @] 21) 
RKY=SFLOATF (KY) -] ,9] 

hte 

Le) 

OF| UMX2N VG 


30 


Bry a) 


an 


} 


Sar 


50) 


Z 


a) Ss (EO 

2 A= OO 00. 

TAUSXTES 

OY 500 THle20n 

U(]) =v 

Taree lL vat, 0 

=o 

XSTGMA()) S020 

i ec) Se 

rue) =1 0 

RXSTOMA(?P)=020 

PAM (a) =0,.0 

xt 3) =O .0 

XSIGMA (3) 210% 

YEXSTAFT 

YSYSTART 

STEPSHTAUS ENO, 

TSHR/24 . 

CPA S02 THl,209 

CALL TERP 

KROSEGAMA thy 

Cé4LL SHIP 

wo ss Wy) 

CALL VDF RIV 

0 550 Jale3 

XE aM (J) STIL AM( YJ) e(CAPVX#COSF (UCT) ) TL AM (J) *CAPVX#SOFNE (UCT) ) XM (YQ) 
+GAMA#EYPF (ARG) HHX#XYSTGOMA(U) ) #STEP 

XM (SV SEXMIC I) (CAPVYFCOSF (UCT) FTLAM(J) *CAPVY®SINE (UCT) ) ®¥MU (YY) 
*GAMAFE XPFE (ARG) F#HY*#XSTGMA(U) )#FSTER 

Th OMe) =x Dares) 

X=X*#CAPVRCOSF (I(T) )HSTEP 

YEY#CAPVRSTNE (UCT) ) #STEP 

he KET A491) EO TN 524 

TF (Y eh Te Pet} IGN TO SPs 

IF (X eB6Te RKXIGN TO G?P4 

JFey ,GT, FKy)GO TN S24 

TEST+STFP 

CONTINUE 

TMIN(L) 2=xMit¢]) 

THU A) Sketi(2) 

SOIDFT= XLAM( 1) FxmMi(e) XE AM (2) Xm (1) 
TLAM())=xXMI)(2)74AD JDE T 

AMUN (}) sexXmiI (2) A/ANIDET 

X¥STGMA(1) =020 

TEAM (2) =e XLAM (A) SAD INE T 

rit ( AP) SX} AM(])) SADINET 

XSTGMA (2) 2000 

TILAM(3)=(°X%tL AM(2)#xXMU(3) —XLAM(3) *#TMI(2) ) AADUDE T 
XMII0 3) = (XLAM (3) *TMU( 1) XL AM( 4) *#XM1 (3) ) AANGDET 
XSIGMA(3)=21.9 

X=XSTART 
=YSTAFT 

Za), 0 

TathikRse4, 

XTFStOen 

71120,0 

717?=9,0 

Z13=tYo0 


Sul 


Si 3 


ra le) 


} 


} 


l 


ae 


awl 


22eC=a=N,0) 
129=0,.0 
70 

NO 5SN3 T2led00 
CALL TERP 
KPESCGAMAtH 
mile) ) tle, 

TES (TV SxTesS 
CALL ANGIE 
“w7(T)=xX1G 
AL (1) =x! T 
wa (J )axk 
CALL SKIP 
wos Ut(T) 
CALL VDERIV 


GRADK (1) =ThLAM(1) #(CAPYIIM#COSF (1101) ) eCAPYV*SINE (t1(T)) ) © XME(1L)*% 


(CAPVUFSTINE (LCT) ) 


* CAPV¥COSFE (UI(T)) ) 


GRADY (J) STL am (2) #(CAPVUFCOSF (UCT) )@CAPVASINE (U(T)) ) © Xmit(2)% 


(CAP VIIESTME (UCT) ) 


¢ CaPyveCOSF (U(T)) ) 


GPab7(T)=STLAM(3) #(CAPVUSCOSF (UCT) ) @CAPVRSINE (H1(T)) ) * XMU(3) # 


(CAPVIEIFSTINE (UCT) ) 


¢ CaPVFCOSF(UCT)) ) 


71127211] +GRADy (1) #GQANYK (1) eSTFPR 
71 2=212eGRADX (IT) #GRANY (1) #STFP 
£13=2713+GRaDX (1) #GRADZ (1) #STEP 
720 a72!)4+GRADY (1) eGPANY(T)#STEP 
72 3=7234+GRADY(T) #GRANZ (1) #STEP 
733=733+4GPAD7(7)#GRAN7Z(1) #STEP 


De) ye) Rs: 


xLaAM( J) eTLahM (J) =(CAPVX#COSF (1I(CT)) *# TL AM( I) ¢CAPVX#STNE (UCT) ) #XMtICU) 
+GAMAHFE XPF (ARG) #HX#XSTGMA(J) y#STEP 

XM (J) SxXri Cd) = (CAPVYFRCASF (UCT) ) FTLAM( J) #CAPVY#SINE (UCI) ) @xXMit( J) 
+GAMA#EXPF (ARG) HHY*#XSIGMA(J) )#STFP 


TE AM( Jy S=xL AM (J) 


K=X*CAPVHCOSF (LI(T))*#STFP 
YEY*CAPVESINE (1(7)) #STEP 
7=Z + FXPF (ARG) #STEP 


T=T*STEP 
XTFSHZXTFSeSTFP 
CONTINUE 


XNOTSCAPVECOSF (LUC 200) ) 
YOOTHSCAPVASINE (U(200)) 


ARG=GBAMAeK 
ZNOTSEXEF (ARG) 
CALL TERP 

Hr (201) =+ 

TES (201) =XTFS 
CALL ANGI F. 

AZ (2N1)exl.6 
Ai (201) =XLT 
wh(2N01)=2XkK 
YIVJFR SXF ANDY 
YOJTFFSYFND=Y 
NFL xsyeF Ox 
NEL YSYor OY 


MSFZ(DELX##C ene y##2+973,75) /1043,638743 
FRROR]=XNIFFaeee + YIIFF oH? 

IF ( ERROR=MSFReMSFHCRITHCRIT#O.00002366 ) 505¢505451]1 
ALL=EZL1L&#2 +71 2HePeKXNNTHRC 


a2 


ALPS711#71°2 ¢2Ze2e71 2e¢xKDOT#YDOT 

BP PBTILCRHHEC ELA CHHC CYDOTHHC 

ME TSALL HACP HA] ote? 

CLS(A2Z2P#XNTEFeALA#YNIFF) JDET 

C2H(ALIL#FHYNIFFeAlLeoe#XxnIFF) JNET 

El=Z711#C1+712%C2 

FP=Z12#C) +72A0#C2 

DEL TAUSCTHEXNOTeC24YNOT 

TAUSTAUEDFEI TAU 

I. SNH Ksle?I0 

WELUCK) SFLYGRANY (K) +F2*GRANY (K) 
GSO 6(K) SUCK) @NELUCK) 


NON + } 
Pare GT. P5600 TO SOT 
aon | Awa) 


SO7 PRINT 9179 
Slo FORMAT (/5X14HNO CONVERGENCE) 
50 TO AS 
SAS CONTINUE 
L=l¢] 
IF (NNeEM.I) CBO TN 544 
1 OA SEP Se PaO L HYG OPO s 1 2 
NPL UNXSPELIIMNXS2, 
Perc RmIMY .UTe 0.01) 460 To 514 
312 TAUOLNSsTaAu 
GA TN S44 
| 544 IF (ZeL T2011) GO TO 545 
EL UMX=SNFILUMX/2.6 
heaping aT. 0.0%) GO TO 6)4 
545 701.022 
546 K11=711%71) + 712%712 «© 713#713 + XDATHEXNOT 
Hi2s7212#711 +¢ 2128722 + 7238#7)3 © XDNTHYNOT 
3271342), ¢ Zaz*Z)2 + 233%7)3 + KOOTHZNOT 
Hees7Zlee712 + 7AA%72P2 + 223#723 + YNONTHYRHOT 
R2327134712 + 723%722 « 733#723 + YDNT#Z0NNT 
A33=713%7)3 ¢ 723*223 + 72334733 ¢ ZOOT#eZNOT 
D=SRLIF(RO2H KISH 4239 4%R23) eRlL2%(R]2#R332R238R13) 
] «RB13# (BR) A#h23.R228813) 
D1 =(B] 348228124823) /D 
Ups (R114#h23=-45138R1 2) /0 
D3e(nblevbh)eeB) 1 #827) sb 
HL=Z114#D] + 712802 + 71 3%#D3 
WPSZ1L2#N] + Z22*hP + Z223#D3 
r3=713#D)] « 7238#NP «+ 7334N3 
B4=2=XNOT#D] ¢ YONOT#H#D?S + ZNOT#N3 
B1G6G=0,0 
NY 515 JTHle2N0 
OFLU(T)SRI#GRADX(T) «+ R2e#GRANY(I) + R3*¥GRADZ (TI) 
SREABSF( DELI(T) ) 
IF( BB eGTe RIG )RIGZBR 
S1S CONTINUE 
RO = DPELUMX/( BIG ¢)]-0*#ARSF(R4)/TAU ) 
09 520 Jsle200 
NELUCT) =PR#DELU(T) 
S?P0 UCIISUCT) «¢ NELU(T) 
IF (L «GT. 25)60 TO S52) 
69 TQ 50} 
5134 CONTINUE 
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letniNeeeat. ff 360 TO) SAa0 
PRINT S2e 

B22 FORMAT (////29X}, AHMINIMUM TIME ROUTE) 
PRINT 5436!) 

S43 FORMAT(S/20XTHNET B =k6e?2) 
PPINT JA 
PRINT 190 ( TES(K) 0 A7(K) sAL(K) OHT(K) eWA(K) 9K21920} 910 ) 
BAmMAS0,01 
ayaa ea | 
GY TO S41 

540 PRINT 542 

S42 FIPMAT(/S///PRX | GHWEIGHTEN TIME ROUTE) 
PRINT 543¢1) 
PAINT JA 
PRINT 3908 ( TES(K) eA7(K) oAL(K) SHT(K) WACK) 9K21920) 010 ) 
GO TO &3 

SP) PRINT 523 

S?3 FARMAT(SXRHMO SOL N) 
69 TN @3 

SPE PRINT SPSaK OY 

Get FORMAT(SX4QHHSHTP TS NOW OUTSTNE THEATER OF OPFRATTONSX3HX =FINeS» 
1 SX3HY =F190e5) 

ma STOP 

FID 
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SURROUTINE TFERP 

DIMENSION COG) 9P (4) 9914) »PX(4) sQY (4) «HP (4) 9 CP (4) 9 SP (4) pHPX (4), 
THPY (4) CPX (4) sCPY (4) 9SPX(4) 9SPY (4) eHN(4) 9CD(4) 9S (4) HS (4) 9CS(4) 9 
259 (4) @HT (464) 9 CT (494) 6ST (404) 

COMMON Xe VoAsReoCCoHoCK eSKeKXeKY oKT oe FOX eFOY e KCON 

COMMOIN/L 1 ZXHT (8448) sCSKk (8448) e SNK (8468) oT 


ee. 
Lie 


COMMON/L2/HX @HyY 
COMMON/IL3/0KX eNKY 
L=XF TXF (T) 

TT=(m INTE (T) oT) #2 6), 
TP] =TT+1,. 

TMl=STT=1, 
TPMSTP]4TM] 

IF (L=KT4#3) 19104 

IF (L)2e92e3 

K4=3 

TM3STT #3, 

Cer STM isiig7s, 
C(2)saTP]#TM3/4, 
C(3)=TOM/B. 

690 TT) 16 

R424 

G=(3e*TT+2.) #TTH9, 
Fae4 HTT 4G 
C(l)seTAieTMjys 16, 
C(2)2=G6*T41/16. 
C(3)seFHeTpl/16, 
C(4)=2TeMeTPlyla, 

COr On 1S 

IF (L=KT)69505 

x42] 

LeskTea] 

Cll)=l. 

GO To 16 
Geo(3.#TT+2,) #TTH9, 
rSeG HTT 4G 

Cll) s=TeumetTmMis1a4, 

TF (LeKTo1)11,10,5 
K4=2 

C (2) =s(G#TM] + (T2MLF)#TP]) 7/16, 
G0 TA 75 
C(2)=GR TM] /1L6- 

IF (L*KT4#2) 3912619 
K4=23 

C3) S(T2M=@F)#TP1/16., 
L=t=-} 

MsXFTXF (X) #2 
MSXFIXF(Y)=2 
XX=(mINTE(X) +X) #2, 0], 
YYS(mINTF(Y)+Y) #2 0], 
API=XX4+], 

XMLEKXel, 

YPJ=yYYel, 

YY HVYVY=1,. 
X2@MEXP1#Xm1 
YOMsyPleym) 
P(]1)=—XM) *X2M/16, 


o5) 


]& 


fel 


ae 


3 
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PA 


P(P)H((3,%XX 42.) *XXH9.) FKMISIA,W 


P(3)2—exx#XX/8e4]1,125eP (2) 
P(4)=2XP]#XOM/ 16.6 
Q(Ll)s—-YMIl#YemMslLA, 


Q(AQ)S0(3e#VV¥ 42.) #YVYe9~.) HYMI1/16. 


W(3) SeV¥VY#VY/A.412.125=0(2) 
(4) 2YPIl#yYemMs1l6, 

PX (4,2 (N-3754%XX=9,.1] 25) #XP] 
PX(])S0.5*%#XKX—-PX (4) 
P¥(2)2(1,l2e>5%x xo], 375) xp) 
PX (3) se S*XKePX (A) 

NY (4) 2(0.375*%YY=9.175) # YP] 
WY (L)SO.5# YY ey (4) 

Oy (2)2(1 ,125¢yyn) , 375) #YP)] 
QV (3) Sen eS*YVHQY (2?) 

lyf) 27 Kol eK4S 

HP (K) 20,0 

CP (K) S020 

SP (K) S020 

HPYX (KSA LN 

HOY (KSC .N 

I SA 

COY (kK) SND 

SPX (K) S09 

SPY (Kis), 0 


KKs( (K#L)#KYON )#KX & M @mKCON 


jaca Let 

RA (Js 0 

CO a Ore 

$9) (J) 20.0 

it) a0 0) 

CS(J) S049 

S$S(J) 20.9 

JJSsJRKXeKK 

DO 23 Ts] 4 

flsJ+JJ 

HT (LleJ) axel (It) 
Cie = 6S (1.1 } 

ST (TeJ) SSNK (TI) 

pa eS T=} o & 

Tes a= ie 

md( IT) SQO(U) FAT (Ted) +H CT) 
CWC Te OS ACT Cle) eC) (1) 
SN (IT) S009) *ST (Led) #Sp(t) 
AS(JTISP CJ) FRT (Ye T) +HS(L) 
CS (Il )y=P (yee) (Us bie 615 
SS(I) =P (J) #ST (Je) +99 (7) 
CONTINUE 

NN 27 T=led 

HP(K) SH (I) ®P(].) +HP(K) 
CP(K) HFCN( 1) FP (J) +cP(K) 

SP (Kj) 25D (1) #P(1)+SP(K) 
HOX(K) SHH (IT) #PX (7) ¢+HPX(K) 
CPX(K)=SCN( LI) #PX(T)+#CPX (K) 
SPX(K) =SH(T) #PxKX (1) +SPx(K) 
HPY(K) SHS (LT) #OY (1) +HPY (K) 
CPY(K)ECS(T) #QY(1) +CPY(K) 
SPY(K)SSS(T)#OY(1I)+SPY(K) 
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aif 


er 


co” 


3() 


=| 


S82 


CONTINUE 

H=0,9 

CK20,0 

Se=0. 0 

HxX¥=0,7 

HY=0.0 

CKXx=0,.0 

CKYS020 

SKX=S0 20 

SKy=0,0 

OO 3] Kei )eKs 

hH=C(K) #HP(K) oH 
CK=C(K) #CP(K) «CK 
SK=C(K)*#SP(K) SK 
HX=C(K) #HPX (CK) HX 
Hy=C(K) #HPY(K) eHy 
CKX=C(K) #CPX(K) +CKX 
CKYESC(K)#CPY(K) ecKy 
SKX=C(K)#SPX(K) +SKX 
SKY=C(K)#SPY(K) +SKY 
CONT TINUE 

KAN=SSORTF (CK&#CK+SK#SK) 
CK=CK/RAD 

SK=SK/RAN 

DK X=CK®SK X=SK#CK xX 
DKY=CKFSKY@SK#CKY 
KE TURN 

END 


oy 


SURROUT Itt eSHTP 
COMMON WeYeAoReCCeHeCKeSK eK XeKVY eK TeSFOX Se FOY eKCON 
COMMON/L CZ HX a HY 
COMMON/LG/SAX ae AY ae BX oe RY oC Xe CY 
SPFEI) PARAMETFIES OF A APY 
FAC] sti+pu. 7) G8HRN79 
VASFACL HC, 49RD ABR PERE XPF (9909347059 30444H) #1 075341054963 
FACL2£2.14949466258/F AC] 
FAC CH=FAC) $0 eNVI47TNHNS9IZIU4SS 
DOVASVHEF AC? 
FAC] SHet ,9AS4HH2PT4 
VE SFAC] #40, 34A431 AP) 9O9BHE YPF (0.025 32469547 364H) #7 96379908524 
FACILHN, 34H 31LC19GAGA/SEF AC] 
FACCH=FAC] #0 .07°53749547 36 
DVESVF#F ACP 
FaC1=H4+ 372997335271) 
HF AC] #354274 360591 HF XP (=0 NRHN 2945246904H) #0 40R4BR73291 Fad 
FAC LS 3 2.58274 36158) J/FACI 
FACC HF AC] Hf ge VHABNAQG4524G9() 
Iranser ACP 
A= (VE OV) Hi eS 
COC=heVH 
tp)As (PVE +P Vid) 0S 
OCsV4e)Vy 
AaX=sDateex 
nYtDAdHy 
mK SOR HH YX 
MVYED3FHY 
CY¥sVUC#hY 
CY¥=DG7 Hy 
RO TUR 
rity 
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SUBROUTINE VNERIV 

COMMON) XeYeNeBReCCebeCK oSK eKX eKY eK TeFOXSFOY eKCON 

COMMON/L3/0KX 8 DKY 

COMMON/LG/4OX ep AY e BX eRYV AEX eg CY 

COMMON/LS/CAPV eCAPVX eCAPVYeCAPVUewW 

REAL MSF eMSFXeMSFYeNUM 

DEL XaXmFOX 

DFEFLYSYoFOY 

MSFS(NEL X##2+NEL Y4#4#24973,.75) /1043,63R8743 

MSF xaDELx/521],8193715 

MSFYSNELY/S2] 8193715 

TE( LR oF Qe 0 )GO TO 1 

COST=H=COSF (w) #CK + SINF (w) #SK 

SINTS=SINF (W)*#CK = COSF (W) #SK 

GO TN @ 
1 COSTSCOSP#CK «+ SINP*#SK 

STINTSSINP#CK = COSP#SK 
2 ACMCCEAM*hA@eCCHCC 

ROONTSSARTF (Be#2% (COSTH#2) sACMC2*#(SINT@&2) = ) 

NMiIMSReAPMNCA/SA 

DFE NOMSBRHCCHCOST/A+ROOT 

VENUM/DENOM 

VGRID=V/8,56604)]6667 

CAPV=EVGRID*MSF 

AVYCKXSAHAK@eCCH#CX 

AMCY2A#AxsCCH#Cy 

FX=(CCHRK+BFCKXMPRHCCHAX/SA)/SA 

FY=(CCPRY*H#CYeRECCHAY/SA) /A 

RPCMA2G@SReRseCCHCCuoAtA 

TERMX= (COSTPFX+BACCHSINT#DKX/S AS (AMCX#SINT#SINT+ReRX*COST#COST + 
1 BPCMAZ#SINT#COST#DKX) /ROOT) #(A*VGRIN/(BRAeMce2) ) 

TERMys(COST#F ysReCCHSINT#DKysAs (AMCY#SINTSeSINT,ReRyeCOSTH#COST 
1 *BPCMA2*SINT#COST#DKY) /ROOT) #(A*VGRID/(R*A2MCO) ) 

VXSVGRID#(2e*#AMCX/SACMC2!] *+RBX/RHRAX/AMTERMX) 

VYEVGRID#(2,%AMCy/A2MC2+By/BeAY/AaeTERMY) 

CAPVXSVX#MSF +VGRID#MSEF X 

CAPVYSVY#MSF 4+ VGRINMMSFY 

VULVGRIDs (BeCCaSINTsAsBPCMA2HSINT#COST/ROOT) /(BaCCaeCnSTsA 4RONT) 

CAPVUSMSF#VU 

KE TURN 

FND 











Bo 





SURROUTTINE ANGLE 


COMMON XeVYeAeReCCeHeCK ep SK eoKX eKY eK TeoF OX eFOYOKCON 


COMMON/I.6/7 XK eXLGeXLT 
FL X=XmF OX 

OFLYSYoFoy 
COSXKSHeDEL XK eCKEDFILY#SK 
SINXKSDFI X#SK MIE lL Y*#cK 
IF (COSX¥R) Col ee 

xKe STONE (90 ee SINK) 

sO TN «A 


XK ZATANE (STNXK/SCOSYXK) #57, 2957795] 


TE (COSXK) 39406 
LF (STAIXK) 5S e4eG 
XKZXKOI]RD, 

59 TO 6 
KXK=XKeJ Rl), 

TE (XK) 7 obi gk 

xK = 340, +XK 


XTSUDFI Xe GASRNTTISSDEL YH, 1 7364818 
YTSVEL X#.) 736GRJARSDELY%29R4BN7T75 


RAN=SOAKRTE (XTHXT © YT#YT) 
XT 059s.) 0 

XL GESSTONFE (90 eeYT) 

BA TO )4 

KIL G2OTANF (YTs XT) #57. 72957795) 
TF (XT) 11014914 

IF CYT) 139120)? 

XLGexlLG © LR. 

G0 oleMmy Jed 

xLG=eyLG@ =| JAN, 


XL TE=e-ATANF (RANIIZ3] 6205) 41142591559 + 99, 


Re TURN 
FANI{) 
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TOENT OTGB 
PROGRAM LFNGTH 


RLOCKS 

PROGRAM LOCAL 

ENTRY POTNTS 
NN00000 NTGR 


| EXTFRNAL SYMBOLS 





: TREDTG 
| FNTRY OTGR 
% FORPTPAN CALL TO THIS PROGRAM===CALL DTGB (A,B,C) 
| t meeWHERF A © ANDRESS OF NTG WORN(LEFT JUSTIFTEN)B1 
: Fs B = ADNRESS TAy B2 
+t C © ADDRESS NEW NTG WORD(LEFT JUSTIFIFO)B3 
+ RETURN WITH NEW OTG IN € 
NTG6A RSSZ ] 
SX6 B83 
SAA T 
SA] A} 
1X] 4&B RGT ADJUST DTG 
| SAa Re 
RJ zXTRFENTG 
| SA} T 
| Sh3 X] 
| LXxX6 l2 
SAG B3 
4 JP NTG6B 
i NATA 0 
| END 
: UNUSF!) STORAGE 21 STATEMENTS 3. SYMBOLS 
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CLOENT > FE TEENS x 


2ROGRAM LENGTH 

AL OCKS 

DRAGRAMS LoOCaL 

-NTRY POTNTS 
onnn00 FIFLNSx 


-XTERNAL SYMROLS 


~¥ DEL T wAR64X RDMST 


FNTRY FIELNSX 


9 

4 PEVISEN BY S, 
ty 

-. 

+ 

oy FORTRAN CALL 
=) mee WHEOF A 
+ 8 = 
be C= 
> 

af Ds 


FIFILMSx &SS7 i 


We SEI FRINGE JUN 1968 MELLONICS 


SELECT FNWF FIELD FROM MASTFERe LINPACK AND FLOAT DATA 


TO THIS PROGRAMewe~ CALL FIELDS (AsReCoD) 

ANNRESS OF 3989 WORD RBIN(CONTAINS FLTOD NATA ON EXIT) 
ANDRESS OF NATE@TIME GROUP WORD 

ANNRESS OF PaRaMETER SELECTION INnEX (0910293) 

WHERE 0 = CDel = CHy@ = U CURR»,3 « y CURR 

ADNRESS OF READ*MASTER ERROR YNDICATOR WORD 

(le FTELD NOT FOUNR OR PARITY ERROR (ph) #1 ON RETURN) 


t SAVE FORTRAN CALL PARAMETERS 
SAP Be DTG 
SX6 Bl 
RX7 X2 
SAK PARAM 
SA7 TABLE +] 
SA3 B3 STORE SELFCTION INDEX 
RXA X3 
SAG COUNTER 
SX7 BS NOT*FOUND INDICATOR 
SA7 PARAMs] 
> PE AM ONF FIFLD FROM MASTER TAPE 
RGN SA2 COUNTER 
SA] NAMTAR«a x? BCD PARAMETER NAME OF FIELD 
RX4 X] 
SAG TABL.F 
Sk] A6 
RJ =XROMST 
7 IINPACK FIELD 
Gad R0 
Sd] PAF 
SAea PARAM 
SB? x2 LOC. UNPACKED FIELD 
S43 T60 1R 
PJ 2zXWAB64x 
tt SHTFT DATA INTO LOWER 68 AITS 
SB] 20 
Sho 3988 
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Srilr T 


B2eBleSHIFT 


CONVERT FIELD TO FLOATING POINT 


$81 
SBe 
SA} 


SA2 
SB3 


TAR| F 


NAMTAB 


SCAT AK 


COUNTER 
PARAM 
ERR 


PAF 


UNUSFD STORAGE 


x1] 

jaa 
COUNTER 
SCATABsey] 
X2 

760)}8 
aXFXDFLT 
FIELDSX 


0 

0 

30/PAF y3n/ERR 
SLTAPES 

3LEEE 
10HCN 
10HCH 
1OHU 
1oKV 
4} 

4] 

36 

36 

0 


CURR 
CURR 


e 
PARAMe] 
l 

Bl 

X] 
FIELNDSX 
1100 


85 STATEMENTS 
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SCALE FACTOR 





TOFNT SEALAND 
PROGOQAM LENGTH 


BLOCKS 

PRIIGK A Mt LOCAL 

FNTRY POINTS 
ONDNVOO SEALAND 

EXTFRAAL SYMHOALS 

LANISF Sh 


FNTRY SEAL AND 
rr FOOTRAN CALI TO THIS PROGRAM#\@= CALL SFALAND(I,JsLAND oKER) 








“WHERE T = apnRPrFSS FO REaL Y INnEX ON FNWC 63X43 GRIN 
J = ADDRESS OF REAL JU INNEX 
¢ LaNn = ANDRESS OF LAND/SFA INOTCATOR (9 FOR SEAe 1} FOR LAN 
ie KER = aNnRESS OF FRROR INNDTcaTOR (OUT OF ROUNDS) 
+ (O IF ERROR ARSENTe | IF ERROR IS PRESENT) 
5 og 
SEAL AND R857 ] 
ay Rl I 
SAo Re J 
IX] X19) 
IX? KX2—9R2 
Sky R1+39 
SH? B2+15 
| X] XleR1] I $39 
x? y¥2e0R2 J $15 
maxX 36 
aX XQ*X} | 
RX? ~xXO0#X? ) 
Ay ] ylex2 I BITS 24.4794 BITS Nao23>% 
SXé A 3 EACH SCALFD S}5 
SxX7 R4 
SAG LSIN SAVE L/S TNDe ANDRESS | 
SA7 FRRYIN SAVE ERROR INDe ADDRESS 
SR) SEA 
SK2 L ANN 
SRA FRROAR 
te =XLANNSFA 
Sk A Xb 0 
mx 0 
yP OUT 
L ANN SX A 1 
MX7 0 
JP OUT 
FRROAO SXé ] : 
SX7 ] 
Ory SA] LSIN | 
SA2 ERRIN 
SA6 X1 
SA7 X2 
JP SEA, AND 
LSIN S67 ] 
ERRIAN KSS7 ] 


F NT) 





DIMENSION OF WAVE FIFLD ARPRAYS=( 22,5 329 12) 11 home dup 
ROUTE OF SHIP REGINS 0 HOURS AFTER TIME SERIES ORTGIN 
FROM LONGIJTIDE = 154,0 AM LATITUDE = 41.0 

TO LONGITUDE = -123.n9 ANID LATITUDE = 3A.0 


LPNS 39 


GFONESTC ROUTE 





DAYS LONGI {A4Tis= WAVE WAVE DPTRECTION 
OF TRAVE| @TIIDE TINE HEIGHT FROM NORTH 
0.9Nn 154.0 41.1 Ple2l ae 
Eee 156.6 42.0 elec 3n4 
1e0n }SS.9 43.9 [Geo 23) 
ib 96 163,8 44.7 )Sac3 25? 
oO Mae ia VRS) ES 2000 27{) 
o Ais ie 171.0 45,4 19.54 19% 
390 dS ie 16678 2645 
3.50 eo 47.) 11244 ao 
4.00 ~1/5.4 PT 9.7} 199 
4.5n -}70,3 a ers | 19,364 208 
a 0n =O e ig Tk 4.09 1? 1 
aia ~160.0 47.4 6.55 23? 
AON -155,0 47 9.48 1S0 
‘ee #“i150.1 G6 .3 11.78 7 
0 -145,5 45.4 17,28 174 
7.O0 -141.2 44.4 13.41 1462 
Kean -~136,8 43.¢ 7.39 153 
etry () -1327,6 41,% 4,86 162 
2.09 =1PR.5 4).c y 94 ?9 
a0) -124,6 4A 5 6.44 344 
9.7n -l123.1 37,8 13.70 314 
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MINIMUM TYME ROUTE 


DET « = N09 


OF 


NAYS LONGI | ATI = wAVE wAvVE DIRECTION 
TRAVEL. “TUDE TUNE HEIGHT FROM NOOTH 
A200 154,0 41.0 21.2) See 
0 4A Der eay G04 14.40 ] 
2o6 161.6 4(ve3 6292 ee 
1244 L662 40.5 10.10 29) 
lie?? 179.2 40.9 14.79 159 
CeA4 Nf 174.3 4). 15264 78) 
2eARA ata I ueel 17-50 280 
3e 34 -1]78.0 44.0 13-88 27s 
3.84 -173,.8 45.1 12649 eRe 
4.3? -169.3 45.9 12-4 275 
4.AN -l64e7 4A, |) 7206 14e 
5e PR “1690.0 4526 6eA9 334 
my. 7h -155.6 4429 9.35 147 
6224 =l5)23 44.3 19.78 157 
be72A -147.3 44.7 13.9) 1464 
T.2Nn -143.4 43.1] 12.97 17) 
7.62 -139.3 4A Dd 9.42 166 
R.16 -)35.1 41,7 5eA2 296 
A464 e139.,8 49.7? 3296 115 
9.19 -)26.8 39.44 1.29 l 
9.6Nn -123.0 38,1 15.20 B 


DAYS 


TRAVE| 


0.09 
0 4H 
aon 

1.45 

1.94 

2242 

2290 

3,39 

3,87 

4.34 

4.84 

5-3? 

5.8) 

629 

6,72 

7.26 

7.74 

R23 

R77) 

9.2n 

9.6,R 


LONGI 
= 1 Ue 


154,0 
lS726 
161.6 
144.0 
Lise 
174.2 
lias 
-178.,4 
-17/4,3 
-149,7 
-)64,.9 
“149.1 
-155.6 
-151.3 
-~147.4 
“143.5 
-139.4 
-135.1 
-139.9 
“126.4 
el22.9 


sFIGHTEN TIMF ROUTE 


-T RS 


| aTI~q 
THRE. 


ie ei 
397.5 
39,8 
4(1),4 
Alice) 
42.9 
44,4 
45,9 
4A,8 
47,9 
44.4 
46 Sie SD 
44,7 
43.8 
aa. 
42.4% 
41.7 
40,7 
39,5 
6 ala 


WA VF 


HEIGHT 


Cleel 
ls. 3? 
5.18 
9.10 
13.90 
Lue Si 
17.61 
13.48 
luis 7 1 
11.14 
5244 
6.41 
9.76 
19.94 
14.09 
12.77 
9.14 
527% 
3,13 
1.1? 
15,90 
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WAVE DNTRECTION 
FROM NORTH 


3)9 
39 
}45 
215 
ag | 
7A) 
2R) 
270 
275 
CAA 
136 
309 
148 
1558 
1465 
17) 
167 
1A7 
109 
359 
183 





GLOSSARY 


TEST 

IMIN JMIN = Minimum indices of wave height 
data subfield 

|, aS Pa Gt = Dimensions of wave data sub- 
field 

DTG (KT) = Array of octal-date time groups 

U(200), DELU(200) = Control arrays 

GRADX, GRADY, GRADZ = Gradient arrays 

AZ, AL, HT, WA = Arrays for track position 
tabulation 

WE (3989), WD(63,63) = Unpacking area of core 

XLAM, XMU, XSIGMA = Components of adjoint vectors 
(ag eee 

Moray, = Rectangular coordinates of ship's 
position relative to subfield 
origin 

ce Ce = Parameters for elliptical polar 
velocity diagram Fig. l 

H, CK, SK = Wave height, cosine and sine 
of wave direction relative to 
OXY axes 

BOX, FOY = Floating point coordinates of 
North Pole relative to subfield 
origin 

a, CSK aaSNK = Wave field arrays 

Th = Time in days from first member 


Of tame serizes 


Wee = Time in days from beginning of 
eEack 

XMEAN, YMEAN = Coordinates of midpoint of sub- 
field 

ROX, ROY = Semi-dimensions of a rectangle 


in the subfield within whech 
track must lie 
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LR 


COSLG1, S@NLG1 


COSsLIl, SiG 


COSLG2, SINLG2 


COSLT2, SINLT2 


ale, Vigo Xo. yo 


S12 


ARC 


XoUARI"M, LoLARYT 


XEND, YEND 


COSP, omer 


GAMA 


DELUMX 
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Index to denote a geodesic track 
(LR=0) or an optimum track 
(ER) 


Cosine and sine of longitude of 
initial point 


Cosine and sine of latitude of 
iad tell Od tie 


Cosine and sine of longitude of 
terminal point 


Cosine and sine of latitude of 
terminal point 


Coordinates of initial and ter- 
Minal points Gi'track celacuve 
to North Pole in grid units 


Straight line distance from Xl, 
Wales CO AL 7 Ye 


Great circle distance from Xl, 
Wel tOie ee 2 


Coordinates of initial point in 
grid units relative to subfield 
erilgan 


Coordinates of terminal point 
in grid units relative to sub- 
field origin 


Cosine and sine of ship's head-— 
ing at any point along great 
circle track 


Control program to define 
inicial Conmael ute) 


Weighting factor used in cal- 
culating optimum track, GAMA=0 
minimum time, GAMA > 0 weighted 
time track 


Maximum allowable change in the 
Gontrol u(t) ae any instance om 
time 


Counters to stop computation if 
Convergence is questionable 


STEP = Variable time steno of integra- 


tion 

XDOT , YDOL. 2DO Lr = Time derivative of state 
variables 

MSF = Map scale factor for stereo- 


Giaphic proj ectvon 


Ald; Al2, A22 = Elements of A matrix of Sec- 
Evon Lil 

i, Lie Zia = Elements of Z matrix of Sec- 

Z23, Geo jabs | o\ 9 ROAE 

DET = Determinant of A matrix of 


Section III 


Cl, 1@2 = Elements of C matrix of Sec- 
muon, 1 i 

Hl eee, DEL EAY = Elements of E matrix of Sec- 
~1on LIL 

NN = Index to denote a minimum time 


track (NN=0) or a weighted 
time track (NN=1) 


Bie able vols, Boze = Elements of B matrix of Sec- 

23.,..Biss trem LIT 

D = Determinant of B matrix 

Die D2, D3 = Elements of D matrix of Sec- 
ie nL(epave JE Jeb 

Bee Bid 3 = Elements of a matrix of Sec- 
caom ol LT 

RR = Coefficient of DELU(T) tomlamie 


maximum change in control at 
any time to DELUMX 


era) = Weighting factors for interpola- 
tion in the time demension be- 
tween K4 ordinates 


L = Index to determine member of 
field time series 


DE = T of Eq. (20) of reference [5] 
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M, N = Indices to pick outa) grid 
point data 


XX, YY = XY of equation (18) of 
reference [4] 


Pig imtoo: PxA) = Pl to P4 of equation (19) of 
reference [4] 


Ol) to OC4) = Elements of pt (y) Matiscisco £ 
equation (18) of reference [4] 


PX, QOY = Partial derivatives of P(1) to 
P44) and O(1) tomo 


15 Dem 9 Os Grime C1 ,6 OOF = Partial derivative of H, cos K, 
SKX, ‘oi Sin K 
DKX, DKY = Partial derivatives of wave 


direction angle K 


Vite Ve = Ship's speed in knots in head 
waves and following waves 


AX eraser AY, BY ,mGy = Partial derivatives of param- 
eters of elliptical polar 
velocity diagram 


VDERIV 


MSFX, MSFY = Partial derivatives of map 
scale factor 


COST, TsiNT = Cosine and sine of the angle 
Oo of Figume | 


VGRID = Ship's speed in knots relative 
to the stereographic grid 


CAPV = Ship's speed in knots relative 
to the earth's surface 


CAPVX, CAPVY, CAPVU = Partial derivatives of CAPV 
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ABSTRACT 


With the advent of the high speed digital computer, many 
problems heretofore considered unsolvable for all practical 
purposes are now well within the reach of the anplied math- 
ematician. One such problem is the routing of a ship through 
a time dependent ocean wave field, from one point on the 
earth's surface to another, so as to minimize a cost func- 
ELON Oretae. form g(x;y,t,uU) - 



















This paper considera a numerical solution to the above 
problem. The technique to be employed is known as the method 
of steepest ascent and is attributed to Arthur E. Bryson and 
Walter F. Denham [1]. Although the computer program as given 
in the Appendix is written specifically for a VC2AP3 class 
vessel operating in a described area of the North Pacific 
Ocean, it can be readily modified to accommodate any type 
vessel operating in the Northern Hemisphere. 
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